\(~\)

The entire mouse embryonic spatial data (containing~ 1.3 million cells) was downloaded from Xenium website. Using Xenium explorer, kidney cells were selected and annotation file (cell and clusters ids) for selected cells were imported into R. From entire mouse seurat data, a smaller subset containing only kidney cells were selected, called “xe”. Moreover, Seurat object called “so” containing developing mouse kidney single cell RNAseq data was imported into R to transfer labels from single cell data to spatial Data.

\(~\)
\(~\)

0.0.1 Preprocessing data. Match gene symbols between the two datasets.

library(Seurat)
library(tidyverse)
library(ggplot2)
library(DT)
library(clustree)

\(~\)
\(~\)

cells_selected <- read.delim('~/Xenium/Mouse/Selection_1_cells_stats.txt' )
colnames(cells_selected) <- cells_selected[2 , ]
cells_selected <- cells_selected[-c(1,2) , ]
dim(cells_selected)
## [1] 48087     4
cells_selected %>% head()
##      Cell ID    Cluster Transcripts Area (µm^2)
## 3 nlndahgm-1 Cluster 21         154         81.1
## 4 cjkdmpeo-1 Cluster 21         170        69.04
## 5 cjdgmgfm-1 Cluster 21          32        18.74
## 6 cjdlgclb-1  Cluster 2          73         33.5
## 7 cjckllcb-1 Cluster 27         118        62.26
## 8 cjffimih-1  Cluster 1          74        46.09
xe <- readRDS('~/Xenium/Mouse/so_mmu_xenium_kidney.rds')

cells_selected <- cells_selected[ match( rownames(xe[[]]) , cells_selected$`Cell ID`) , ]

xe <- AddMetaData(object = xe , metadata = cells_selected$Cluster, col.name = 'Clusters')

obs <- xe[[]]

Idents(xe)  <- "Clusters"
Idents(xe) %>% unique()
##  [1] Cluster 5  Cluster 2  Cluster 21 Cluster 27 Cluster 39 Cluster 4 
##  [7] Cluster 1  Cluster 32 Cluster 47 Cluster 41 Cluster 66 Cluster 34
## [13] Cluster 14 Cluster 8  Cluster 56 Cluster 12 Cluster 37 Cluster 17
## [19] Cluster 48 Cluster 59 Cluster 40 Cluster 65 Cluster 15 Cluster 31
## [25] Cluster 6  Cluster 50 Cluster 43 Cluster 58 Cluster 3  Cluster 23
## [31] Cluster 20 Cluster 44 Cluster 52 Cluster 7  Cluster 36 Cluster 60
## [37] Cluster 38 Cluster 29 Cluster 19 Cluster 53 Cluster 10 Cluster 18
## [43] Cluster 28 Cluster 11 Cluster 30 Cluster 13 Cluster 55 Cluster 62
## [49] Cluster 16
## 49 Levels: Cluster 5 Cluster 2 Cluster 21 Cluster 27 Cluster 39 ... Cluster 16
xe
## An object of class Seurat 
## 541 features across 48087 samples within 4 assays 
## Active assay: Xenium (379 features, 0 variable features)
##  1 layer present: counts
##  3 other assays present: BlankCodeword, ControlCodeword, ControlProbe
##  1 spatial field of view present: fov
# Selecting cells having more than 10 molecules
xe <- subset(x = xe , subset = nCount_Xenium > 10)
xe
## An object of class Seurat 
## 541 features across 48033 samples within 4 assays 
## Active assay: Xenium (379 features, 0 variable features)
##  1 layer present: counts
##  3 other assays present: BlankCodeword, ControlCodeword, ControlProbe
##  1 spatial field of view present: fov

\(~\)

so <- readRDS("/mnt/market/anclab-rstudio-server/home/mpir0002/Project1/E18_01.2023/E18_Final_for_GRN.rds")

so
## An object of class Seurat 
## 32994 features across 4951 samples within 2 assays 
## Active assay: SCT (16492 features, 0 variable features)
##  2 layers present: counts, data
##  1 other assay present: RNA
##  2 dimensional reductions calculated: umap, pca

\(~\)

DimPlot(so, reduction = "umap", label = T , dims = 1:2 , label.size = 6) + ggtitle("Original single cell data")

\(~\)

DefaultAssay(so) <- 'RNA'

so[["SCT"]] <- NULL

so <- so[which(rownames(so) %in% rownames(xe)) , ]
so
## An object of class Seurat 
## 295 features across 4951 samples within 1 assay 
## Active assay: RNA (295 features, 0 variable features)
##  2 layers present: counts, data
##  2 dimensional reductions calculated: umap, pca
n <- rownames(xe)
n <- n[n %in% rownames(so)]
xe <- subset(xe, features = n)
xe
## An object of class Seurat 
## 295 features across 48033 samples within 1 assay 
## Active assay: Xenium (295 features, 0 variable features)
##  1 layer present: counts
##  1 spatial field of view present: fov
#save.image("1_preprocessing.RData")

\(~\)
\(~\)

so <- SCTransform(object = so , assay = 'RNA')
so <- FindVariableFeatures(so)
so <- RunPCA(so)
so <- RunUMAP(so, dims = 1:30, reduction = "pca", assay = "SCT" , return.model = T)

xe <- SCTransform(object = xe , assay = 'Xenium')
xe <- FindVariableFeatures(xe)
xe <- RunPCA(xe)
xe <- RunUMAP(xe, dims = 1:30, reduction = "pca", assay = "SCT" , return.model = T)

#save.image("2_preprocessing.RData")

\(~\)

load("2_preprocessing.RData")

\(~\)

0.0.1.1 Generating UMAPs using common genes between the two datasets (~ 295 genes)

so      
## An object of class Seurat 
## 584 features across 4951 samples within 2 assays 
## Active assay: SCT (289 features, 289 variable features)
##  3 layers present: counts, data, scale.data
##  1 other assay present: RNA
##  2 dimensional reductions calculated: umap, pca
xe       
## An object of class Seurat 
## 590 features across 48033 samples within 2 assays 
## Active assay: SCT (295 features, 295 variable features)
##  3 layers present: counts, data, scale.data
##  1 other assay present: Xenium
##  2 dimensional reductions calculated: pca, umap
##  1 spatial field of view present: fov

\(~\)

DimPlot(so, reduction = "umap", label = T , dims = 1:2 , label.size = 5) + ggtitle("single cell data") + theme(legend.position = "none")

DimPlot(xe, reduction = "umap", label = T , dims = 1:2 , label.size = 5) + ggtitle("Spatial Data") + theme(legend.position = "none")

\(~\)
\(~\)

0.0.2 Label annotation with data integration

\(~\)

load("Integrated_data1.RData")
intergration_features <- SelectIntegrationFeatures(object.list = c("xenium"= xe , "so" = so), nfeatures = 2000)

datasets_list <- PrepSCTIntegration(object.list = c("xenium"= xe , "so" = so), anchor.features = intergration_features)

integration_anchors <- FindIntegrationAnchors(object.list = datasets_list , dims = 1:30)

integrated_data <- IntegrateData(anchorset = integration_anchors, dims = 1:30)

#save.image("Integrated_data1.RData")
integrated_data
## An object of class Seurat 
## 1174 features across 52984 samples within 4 assays 
## Active assay: integrated (289 features, 289 variable features)
##  1 layer present: data
##  3 other assays present: Xenium, SCT, RNA
##  1 spatial field of view present: fov
integrated_data[["integrated"]]@data[1:20,1:10]
## 20 x 10 sparse Matrix of class "dgCMatrix"
##                                                                            
## Calb1   .        0.6931472 .         1.3862944 .         1.386294 .        
## Ctss    .        0.6931472 .         .         .         .        .        
## Ptn     2.564949 2.1972246 3.1354942 2.3978953 2.4849066 1.386294 1.0986123
## Podxl   1.386294 2.6390573 .         .         .         .        2.5649494
## Aadat   .        .         .         .         .         .        .        
## Lrp2    .        .         .         .         .         .        .        
## Pf4     .        .         .         .         .         .        .        
## Stfa2l1 .        .         .         .         .         .        .        
## Stmn1   1.386294 1.3862944 2.3978953 0.6931472 1.6094379 1.791759 2.0794415
## Fxyd6   .        1.0986123 0.6931472 1.3862944 0.6931472 1.098612 .        
## Gatm    .        0.6931472 .         .         .         .        .        
## Nts     .        .         .         .         .         .        .        
## Upk3b   .        .         .         .         .         .        .        
## Epcam   .        .         .         0.6931472 .         .        .        
## Cfh     1.386294 1.0986123 2.8903718 .         2.1972246 .        0.6931472
## Plvap   .        1.7917595 .         .         .         1.098612 1.7917595
## Lum     .        .         .         .         .         .        .        
## Cryab   .        0.6931472 .         .         .         .        0.6931472
## Tagln   .        .         .         .         .         .        .        
## F13a1   .        .         .         .         .         .        .        
##                                     
## Calb1   1.0986123 .        .        
## Ctss    .         .        .        
## Ptn     2.5649494 1.098612 1.3862944
## Podxl   0.6931472 2.197225 .        
## Aadat   .         .        .        
## Lrp2    .         .        .        
## Pf4     .         .        .        
## Stfa2l1 .         .        .        
## Stmn1   2.4849066 1.945910 2.9957323
## Fxyd6   1.3862944 .        2.4849066
## Gatm    .         .        .        
## Nts     .         .        .        
## Upk3b   .         .        .        
## Epcam   1.0986123 .        0.6931472
## Cfh     1.7917595 1.098612 .        
## Plvap   .         1.609438 .        
## Lum     .         .        .        
## Cryab   .         .        .        
## Tagln   .         .        0.6931472
## F13a1   .         .        .

\(~\)

integrated_data <- ScaleData(integrated_data)
integrated_data <- RunPCA(integrated_data)
integrated_data <- RunUMAP(object = integrated_data , dims = 1:30)

#save.image("Integrated_data2.RData")
load("Integrated_data2.RData")
integrated_data1 <- integrated_data

\(~\)

integrated_data[[]]$Clusters %>% unique()
##  [1] "Cluster 5"                         "Cluster 2"                        
##  [3] "Cluster 21"                        "Cluster 27"                       
##  [5] "Cluster 39"                        "Cluster 4"                        
##  [7] "Cluster 1"                         "Cluster 32"                       
##  [9] "Cluster 47"                        "Cluster 41"                       
## [11] "Cluster 66"                        "Cluster 34"                       
## [13] "Cluster 14"                        "Cluster 8"                        
## [15] "Cluster 56"                        "Cluster 12"                       
## [17] "Cluster 37"                        "Cluster 17"                       
## [19] "Cluster 48"                        "Cluster 59"                       
## [21] "Cluster 40"                        "Cluster 65"                       
## [23] "Cluster 15"                        "Cluster 31"                       
## [25] "Cluster 6"                         "Cluster 50"                       
## [27] "Cluster 43"                        "Cluster 58"                       
## [29] "Cluster 3"                         "Cluster 23"                       
## [31] "Cluster 20"                        "Cluster 44"                       
## [33] "Cluster 52"                        "Cluster 7"                        
## [35] "Cluster 36"                        "Cluster 60"                       
## [37] "Cluster 38"                        "Cluster 29"                       
## [39] "Cluster 19"                        "Cluster 53"                       
## [41] "Cluster 10"                        "Cluster 18"                       
## [43] "Cluster 28"                        "Cluster 11"                       
## [45] "Cluster 30"                        "Cluster 13"                       
## [47] "Cluster 55"                        "Cluster 62"                       
## [49] "Cluster 16"                        "Cortical_stroma"                  
## [51] "Nephron_progenitor"                "Renal_vesicle"                    
## [53] "Committing_nephron_progenitor"     "Ureteric_tip"                     
## [55] "Medullary_stroma"                  "Proximal_tubule"                  
## [57] "Early_proximal_tubule"             "Collecting_duct_associated_stroma"
## [59] "Podocyte"                          "Sshaped_body"                     
## [61] "Immune_cells"                      "Vasculature"                      
## [63] "Ureteric_stroma"                   "Stromal_nephron_progenitor"       
## [65] "Loop_of_henle"
DimPlot(integrated_data, reduction = "umap", label = T , label.size = 5 ) + theme(legend.position = "none")

\(~\)
\(~\)

0.0.3 Reclustering integrated data

\(~\)

load("Integrated_data3.RData")
integrated_data2 <- integrated_data
integrated_data <- FindNeighbors(integrated_data, return.neighbor = FALSE, compute.SNN = T,
                     prune.SNN = 1/15, nn.method = "annoy", n.trees = 200, annoy.metric = "euclidean", 
                     dims = 1:30 , reduction = "pca" , verbose = F , do.plot = F )

# Cluster Determination
integrated_data = FindClusters(integrated_data, modularity.fxn = 1 ,
                   initial.membership = NULL, node.sizes = NULL,
                   resolution = 1 , algorithm = 1,
                   n.start = 100, n.iter = 100, random.seed = 0,
                   group.singletons = T , temp.file.location = NULL,
                   edge.file.name = NULL, verbose = F)       

#save.image("Integrated_data3.RData")

\(~\)

DimPlot(integrated_data , label = T,  label.size = 6) + theme(legend.position = "none")

\(~\)
\(~\)

0.0.4 Label transfer through MapQuery method

load("Query_mapping.RData")

\(~\)

#https://satijalab.org/seurat/articles/integration_mapping
#xe <- RunUMAP(xe, dims = 1:30, reduction = "pca", assay = "SCT" , return.model=TRUE)
#so <- RunUMAP(so, dims = 1:30, reduction = "pca", assay = "SCT" , return.model=TRUE)

anchors1 <- FindTransferAnchors(reference = so , query = xe , verbose = T, reference.assay = "SCT", query.assay = "SCT" , dims = 1:30 , reference.reduction = "pca")

predictions <- MapQuery(anchorset = anchors1, reference = so, query = xe, refdata = list(celltype = "Clusters"), reference.reduction = "pca", reduction.model = "umap")

Idents(predictions) <- "predicted.celltype"

#save.image("Query_mapping.RData")

\(~\)

predictions$predicted.celltype.score %>% head(20)
## aangeepl-1 aangfoon-1 aanggbam-1 aanggimp-1 aangjbii-1 aangkabe-1 aangkncb-1 
##  0.5151041  1.0000000  0.4727102  0.4452291  0.4580433  0.4128006  1.0000000 
## aangldgn-1 aangmeek-1 aangmohp-1 aangndif-1 aangnejj-1 aangnpca-1 aangoagl-1 
##  0.7106872  1.0000000  0.6165119  0.3698561  0.7847812  0.7139676  0.8927964 
## aangobco-1 aangofkj-1 aangohbb-1 aangpcig-1 aanhbbog-1 aanhbebj-1 
##  0.7026645  0.9807297  0.6422881  0.9266697  0.8314558  0.9645079

\(~\)

DimPlot(predictions , label = T , label.size = 5)

\(~\)
\(~\)

0.0.5 Plots

0.0.5.1 Label annotation with data integration

0.0.5.1.1 DimPlot
DimPlot(integrated_data1, reduction = "umap", label = T , label.size = 6) + theme(legend.position = "none")

0.0.5.1.2 ImageDimPlot
ImageDimPlot(integrated_data1, fov = "fov", cols = "polychrome", axes = TRUE)

0.0.5.1.3 clustree
#clustree(integrated_data1 , node_text_size = 4)

\(~\)
\(~\)

0.0.5.2 Reclustering integrated data

0.0.5.2.1 DimPlot
DimPlot(integrated_data2, reduction = "umap", label = T , label.size = 6)

0.0.5.2.2 ImageDimPlot
ImageDimPlot(integrated_data2, fov = "fov", cols = "polychrome", axes = TRUE)

0.0.5.2.3 clustree
#clustree(integrated_data2 , node_text_size = 4)

\(~\)
\(~\)

0.0.5.3 Label transfer through MapQuery method

0.0.5.3.1 DimPlot
DimPlot(predictions , label = T , label.size = 5)

0.0.5.3.2 ImageDimPlot
ImageDimPlot(predictions , fov = "fov", cols = "polychrome", axes = TRUE)

0.0.5.3.3 clustree
#clustree(predictions , node_text_size = 4)

\(~\)
\(~\)